function post_plot
%% ☆☆POST-PROCESSING~~☆☆ ☆☆~START！☆☆
if ~exist('plot','dir')
    mkdir('plot');
end
load result\T_up_596_600.mat T_up
load result\T_down_596_600.mat T_down
load temp\str_data.mat str_data
load temp\T.mat T
load temp\W.mat W
load temp\stru_grid.mat stru_grid
load temp\P_aero_up.mat P_aero_up
load temp\P_aero_down.mat P_aero_down
load temp\Q_up.mat Q_up
load temp\Q_down.mat Q_down
load temp\D_T.mat D_T
stru_grid_out=stru_grid((stru_grid(:,1)>400000)|((200000<stru_grid(:,1))&(stru_grid(:,1)<300000)),:);
W=W(:,end); W=reshape(W,6,size(W,1)/6)';
stru_grid=str_data.stru_grid; aerobox=str_data.aerobox;
stru_up_rigid=str_data.stru_up_rigid; stru_down_rigid=str_data.stru_down_rigid;
stru_quad_up=str_data.stru_quad_up; stru_quad_down=str_data.stru_quad_down;
stru_quad_up_plus=str_data.stru_quad_up_plus; stru_quad_down_plus=str_data.stru_quad_down_plus;
[aero_up_pre,aero_down_pre,~,~,~,~]=interpolate_point(1000,aerobox,stru_quad_up);
Npath='E:\nastran\bin\nast20180.exe';
dt=4;
t_end=600;
%% 表面温度场
f=fopen('plot\T_surf.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)),', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n' ...
    ' VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:length(T_up)
    fprintf(f,'%e ',T_up(i));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:length(T_down)
    fprintf(f,'%e ',T_down(i));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

for i=1:length(T_up)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_up(i),2),aerobox(aerobox(:,1)==stru_quad_up(i),6),aerobox(aerobox(:,1)==stru_quad_up(i),10),aerobox(aerobox(:,1)==stru_quad_up(i),14));
end
for i=1:length(T_down)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_down(i),2),aerobox(aerobox(:,1)==stru_quad_down(i),6),aerobox(aerobox(:,1)==stru_quad_down(i),10),aerobox(aerobox(:,1)==stru_quad_down(i),14));
end
fclose(f);
%% 结构温度场
[T_up,T_down]=surface_T(T,stru_up_rigid,stru_down_rigid,aero_up_pre,aero_down_pre);
f=fopen('plot\T_stru.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)),', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n' ...
    ' VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:length(T_up)
    fprintf(f,'%e ',T_up(i));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:length(T_down)
    fprintf(f,'%e ',T_down(i));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

for i=1:length(T_up)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_up(i),2),aerobox(aerobox(:,1)==stru_quad_up(i),6),aerobox(aerobox(:,1)==stru_quad_up(i),10),aerobox(aerobox(:,1)==stru_quad_up(i),14));
end
for i=1:length(T_down)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_down(i),2),aerobox(aerobox(:,1)==stru_quad_down(i),6),aerobox(aerobox(:,1)==stru_quad_down(i),10),aerobox(aerobox(:,1)==stru_quad_down(i),14));
end
fclose(f);
%% 热流
f=fopen('plot\Q.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)),', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n' ...
    ' VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:length(Q_up)
    fprintf(f,'%e ',Q_up(i));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:length(Q_down)
    fprintf(f,'%e ',Q_down(i));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

for i=1:length(Q_up)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_up(i),2),aerobox(aerobox(:,1)==stru_quad_up(i),6),aerobox(aerobox(:,1)==stru_quad_up(i),10),aerobox(aerobox(:,1)==stru_quad_up(i),14));
end
for i=1:length(Q_down)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_down(i),2),aerobox(aerobox(:,1)==stru_quad_down(i),6),aerobox(aerobox(:,1)==stru_quad_down(i),10),aerobox(aerobox(:,1)==stru_quad_down(i),14));
end
fclose(f);
%% 结构压强场
f=fopen('plot\P.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid_out,1)),', Elements=',num2str(size(aerobox,1)),', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n' ...
    ' VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);

l=1;
for i=1:size(stru_grid_out,1)
    fprintf(f,'%e ',stru_grid_out(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid_out,1)
    fprintf(f,'%e ',stru_grid_out(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid_out,1)
    fprintf(f,'%e ',stru_grid_out(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:length(P_aero_up)
    fprintf(f,'%e ',P_aero_up(i));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:length(P_aero_down)
    fprintf(f,'%e ',P_aero_down(i));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

for i=1:length(P_aero_up)
    fprintf(f,'%d %d %d %d\n',find(stru_grid_out(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),2)+200000),...
        find(stru_grid_out(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),6)+200000),...
        find(stru_grid_out(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),10)+200000),...
        find(stru_grid_out(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),14)+200000));
end
for i=1:length(P_aero_down)
    fprintf(f,'%d %d %d %d\n',find(stru_grid_out(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),2)+400000),...
        find(stru_grid_out(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),6)+400000),...
        find(stru_grid_out(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),10)+400000),...
        find(stru_grid_out(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),14)+400000));
end
fclose(f);
%% 力致节点位移场(z方向）
f=fopen('plot\W_p.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Node UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)),', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n' ...
    ' DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
WW=W;
l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',WW(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

for i=1:length(T_up)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_up(i),2),aerobox(aerobox(:,1)==stru_quad_up(i),6),...
        aerobox(aerobox(:,1)==stru_quad_up(i),10),aerobox(aerobox(:,1)==stru_quad_up(i),14));
end
for i=1:length(T_down)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_down(i),2),aerobox(aerobox(:,1)==stru_quad_down(i),6),...
        aerobox(aerobox(:,1)==stru_quad_down(i),10),aerobox(aerobox(:,1)==stru_quad_down(i),14));
end
fclose(f);
%% 热致节点位移场(z方向）
f=fopen('plot\W_T.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Node UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)),', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n' ...
    ' DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
WW=D_T(:,2:7);
l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',WW(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

for i=1:length(T_up)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_up(i),2),aerobox(aerobox(:,1)==stru_quad_up(i),6),...
        aerobox(aerobox(:,1)==stru_quad_up(i),10),aerobox(aerobox(:,1)==stru_quad_up(i),14));
end
for i=1:length(T_down)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_down(i),2),aerobox(aerobox(:,1)==stru_quad_down(i),6),...
        aerobox(aerobox(:,1)==stru_quad_down(i),10),aerobox(aerobox(:,1)==stru_quad_down(i),14));
end
fclose(f);
%% 总节点位移场(z方向）
f=fopen('plot\W_pT.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Node UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)),', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n' ...
    ' DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
WW=W+D_T(:,2:7);
l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',WW(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

for i=1:length(T_up)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_up(i),2),aerobox(aerobox(:,1)==stru_quad_up(i),6),...
        aerobox(aerobox(:,1)==stru_quad_up(i),10),aerobox(aerobox(:,1)==stru_quad_up(i),14));
end
for i=1:length(T_down)
    fprintf(f,'%d %d %d %d\n',aerobox(aerobox(:,1)==stru_quad_down(i),2),aerobox(aerobox(:,1)==stru_quad_down(i),6),...
        aerobox(aerobox(:,1)==stru_quad_down(i),10),aerobox(aerobox(:,1)==stru_quad_down(i),14));
end
fclose(f);
%% 力致应力
f=fopen(['result\101_',num2str(t_end),'_',num2str(t_end+dt),'.pch']);
z=1;
stress=zeros(0,2);
while true
    tline1=fgetl(f);
    if length(strfind(tline1,'HEXA'))==1
        while true
            tline1=fgetl(f);
            if tline1(1)=='$'
                break
            end
            stress(z,1)=sscanf(tline1(1:16),'%d');
            for j=1:3
                tline1=fgetl(f);
            end
            stress(z,2)=sscanf(tline1(60:72),'%f');
            z=z+1;
            for k=1:60
                fgetl(f);
            end
        end
        break
    end
end
stress(:,2)=stress(:,2)/1000;
f=fopen('plot\stress_p.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)+size(stru_quad_up_plus,1)+size(stru_quad_down_plus,1)-2),...
    ', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_up)
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_down)
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_up_plus(:,1))
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_down_plus(:,1))
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

p=0;
for i=1:length(P_aero_up)
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),2)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),6)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),10)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),14)));
    p=p+1;
end
for i=1:length(P_aero_down)
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),2)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),6)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),10)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),14)));
    p=p+1;
end
for i=1:size(stru_quad_up_plus,1)
    if stru_quad_up_plus(i,1)==-1
        break
    end
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_up_plus(i,3)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,4)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,5)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,6)));
    p=p+1;
end
for i=1:size(stru_quad_down_plus,1)
    if stru_quad_down_plus(i,1)==-1
        break
    end
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_down_plus(i,3)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,4)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,5)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,6)));
    p=p+1;
end
fclose('all');
stress_p=stress;
save result\stress_p.mat stress_p
%% 力致应力ini
f=fopen(['result\101_0_',num2str(dt),'.pch']);
z=1;
stress=zeros(0,2);
while true
    tline1=fgetl(f);
    if length(strfind(tline1,'HEXA'))==1
        while true
            tline1=fgetl(f);
            if tline1(1)=='$'
                break
            end
            stress(z,1)=sscanf(tline1(1:16),'%d');
            for j=1:3
                tline1=fgetl(f);
            end
            stress(z,2)=sscanf(tline1(60:72),'%f');
            z=z+1;
            for k=1:60
                fgetl(f);
            end
        end
        break
    end
end
stress(:,2)=stress(:,2)/1000;
f=fopen('plot\stress_p_ini.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)+size(stru_quad_up_plus,1)+size(stru_quad_down_plus,1)-2),...
    ', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_up)
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_down)
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_up_plus(:,1))
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_down_plus(:,1))
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

p=0;
for i=1:length(P_aero_up)
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),2)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),6)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),10)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),14)));
    p=p+1;
end
for i=1:length(P_aero_down)
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),2)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),6)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),10)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),14)));
    p=p+1;
end
for i=1:size(stru_quad_up_plus,1)
    if stru_quad_up_plus(i,1)==-1
        break
    end
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_up_plus(i,3)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,4)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,5)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,6)));
    p=p+1;
end
for i=1:size(stru_quad_down_plus,1)
    if stru_quad_down_plus(i,1)==-1
        break
    end
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_down_plus(i,3)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,4)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,5)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,6)));
    p=p+1;
end
fclose('all');
stress_p_ini=stress;
save result\stress_p_ini.mat stress_p_ini
%% 被忽略的热应力
fold=fopen('153.bdf','r');
fnew=fopen('153+.bdf','w');
while feof(fold)==0
    entry=fgetl(fold);
    if isempty(entry)
        fprintf(fnew,'\n');
        continue
    end
    for i=1:min(5,length(entry))
        if entry(i)=='$'
            entry(i)=[];
            fprintf(fnew,'%s\n',entry);
            fgetl(fold)
            fgetl(fold)
            entry=fgetl(fold);
        end
    end
    fprintf(fnew,'%s\n',entry);
end
fclose('all');
!del 153.bdf
!rename 153+.bdf 153.bdf
copyfile(['result\159_',num2str(t_end),'_',num2str(t_end+dt),'.pch'],'159.pch');
run_153(Npath)
!del 153.f04;
!del 153.f06;
!del 153.log;
!del 153.MASTER;
!del 153.DBALL;
!del 153.IFPDAT;
!del 159.pch;
f=fopen('153.pch');
z=1;
stress=zeros(0,2);
while true
    tline1=fgetl(f);
    if length(strfind(tline1,'$LOAD FACTOR =   1.0000000E+00'))==1
        while true
            tline1=fgetl(f);
            if tline1(1)=='$'
                break
            end
            stress(z,1)=sscanf(tline1(1:16),'%d');
            for j=1:3
                tline1=fgetl(f);
            end
            stress(z,2)=sscanf(tline1(41:58),'%f');
            z=z+1;
            for k=1:45
                fgetl(f);
            end
        end
        break
    end
end
fclose(f);
!del 153.pch;
stress(:,2)=stress(:,2)/1000;
f=fopen('plot\stress_T_ignore.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)+size(stru_quad_up_plus,1)+size(stru_quad_down_plus,1)-2),...
    ', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_up)
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_down)
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_up_plus(:,1))
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_down_plus(:,1))
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

p=0;
for i=1:length(P_aero_up)
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),2)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),6)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),10)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),14)));
    p=p+1;
end
for i=1:length(P_aero_down)
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),2)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),6)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),10)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),14)));
    p=p+1;
end
for i=1:size(stru_quad_up_plus,1)
    if stru_quad_up_plus(i,1)==-1
        break
    end
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_up_plus(i,3)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,4)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,5)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,6)));
    p=p+1;
end
for i=1:size(stru_quad_down_plus,1)
    if stru_quad_down_plus(i,1)==-1
        break
    end
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_down_plus(i,3)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,4)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,5)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,6)));
    p=p+1;
end
fclose(f);
stress_T_ignore=stress;
save result\stress_T_ignore.mat stress_T_ignore
%% 热致应力
copyfile('result\159_500_504.pch','159.pch');
fold=fopen('159.pch','r');
fnew=fopen('159+.pch','w');
while feof(fold)==0
    entry=fgetl(fold);
    if length(entry)<2
        continue
    end
    if length(strfind(entry,'TEMP*                  1'))==1
        entry(44:56)=' 3.000000E+02';
    end
    fprintf(fnew,'%s\n',entry);
end
fclose('all');
!del 159.pch
!rename 159+.pch 159.pch
run_153(Npath)
!del 153.f04;
!del 153.f06;
!del 153.log;
!del 153.MASTER;
!del 153.DBALL;
!del 153.IFPDAT;

f=fopen('153.pch');
z=1;
stress=zeros(0,2);
while true
    tline1=fgetl(f);
    if length(strfind(tline1,'$LOAD FACTOR =   1.0000000E+00'))==1
        while true
            tline1=fgetl(f);
            if tline1(1)=='$'
                break
            end
            stress(z,1)=sscanf(tline1(1:16),'%d');
            for j=1:3
                tline1=fgetl(f);
            end
            stress(z,2)=sscanf(tline1(41:58),'%f');
            z=z+1;
            for k=1:45
                fgetl(f);
            end
        end
        break
    end
end
stress(:,2)=stress(:,2)/1000;
f=fopen('plot\stress_T.plt','w');
fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
    ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)+size(stru_quad_up_plus,1)+size(stru_quad_down_plus,1)-2),...
    ', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,2));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,3));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stru_grid,1)
    fprintf(f,'%e ',stru_grid(i,4));
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

l=1;
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_up)
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_down)
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_up_plus(:,1))
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
for i=1:size(stress,1)
    if ismember(stress(i,1)-200000,stru_quad_down_plus(:,1))
        fprintf(f,'%e ',stress(i,2));
    else
        continue
    end
    if l==5
        fprintf(f,'\n');
        l=1;
    else
        l=l+1;
    end
end
fprintf(f,'\n');

p=0;
for i=1:length(P_aero_up)
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),2)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),6)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),10)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),14)));
    p=p+1;
end
for i=1:length(P_aero_down)
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),2)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),6)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),10)),...
        find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),14)));
    p=p+1;
end
for i=1:size(stru_quad_up_plus,1)
    if stru_quad_up_plus(i,1)==-1
        break
    end
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_up_plus(i,3)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,4)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,5)),...
        find(stru_grid(:,1)==stru_quad_up_plus(i,6)));
    p=p+1;
end
for i=1:size(stru_quad_down_plus,1)
    if stru_quad_down_plus(i,1)==-1
        break
    end
    fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_down_plus(i,3)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,4)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,5)),...
        find(stru_grid(:,1)==stru_quad_down_plus(i,6)));
    p=p+1;
end
fclose('all');
!del 159.pch;
!del 153.pch;
stress_T=stress;
save result\stress_T.mat stress_T
%% 总应力（弃案）
% copyfile('result\159_400_404.pch','159.pch');
% fold=fopen('159.pch','r');
% fnew=fopen('159+.pch','w');
% while feof(fold)==0
%     entry=fgetl(fold);
%     if length(entry)<2
%         continue
%     end
%     if length(strfind(entry,'TEMP*                  1'))==1
%         entry(44:56)=' 3.000000E+02';
%     end
%     fprintf(fnew,'%s\n',entry);
% end
% fclose('all');
% !del 159.pch
% !rename 159+.pch 159.pch
% 
% WW=W+D_T(:,2:7);
% 
% f=fopen('spc2.bdf','w');
% for i=1:size(stru_grid,1)
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),1,WW(i,1));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),2,WW(i,2));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),3,WW(i,3));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),4,WW(i,4));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),5,WW(i,5));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),6,WW(i,6));
% end
% fclose(f);
% 
% eval(['!',Npath,' 153_post.bdf']);
% check_f06('153_post.f06');
% delete('153_post.f04');
% delete('153_post.f06');
% delete('153_post.log');
% delete('153_post.DBALL');
% delete('153_post.MASTER');
% delete('153_post.IFPDAT');
% 
% f=fopen('153_post.pch');
% z=1;
% stress=zeros(0,2);
% while true
%     tline1=fgetl(f);
%     if length(strfind(tline1,'$LOAD FACTOR =   1.0000000E+00'))==1
%         while true
%             tline1=fgetl(f);
%             if tline1(1)=='$'
%                 break
%             end
%             stress(z,1)=sscanf(tline1(1:16),'%d');
%             for j=1:3
%                 tline1=fgetl(f);
%             end
%             stress(z,2)=sscanf(tline1(41:58),'%f');
%             z=z+1;
%             for k=1:45
%                 fgetl(f);
%             end
%         end
%         break
%     end
% end
% !del 153_post.pch;
% 
% stress(:,2)=stress(:,2)/1000;
% f=fopen('plot\stress_pT.plt','w');
% fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
%     ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)+size(stru_quad_up_plus,1)+size(stru_quad_down_plus,1)-2),...
%     ', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
% l=1;
% for i=1:size(stru_grid,1)
%     fprintf(f,'%e ',stru_grid(i,2));
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% l=1;
% for i=1:size(stru_grid,1)
%     fprintf(f,'%e ',stru_grid(i,3));
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% l=1;
% for i=1:size(stru_grid,1)
%     fprintf(f,'%e ',stru_grid(i,4));
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% l=1;
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_up)
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_down)
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_up_plus(:,1))
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_down_plus(:,1))
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% p=0;
% for i=1:length(P_aero_up)
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),2)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),6)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),10)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),14)));
%     p=p+1;
% end
% for i=1:length(P_aero_down)
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),2)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),6)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),10)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),14)));
%     p=p+1;
% end
% for i=1:size(stru_quad_up_plus,1)
%     if stru_quad_up_plus(i,1)==-1
%         break
%     end
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_up_plus(i,3)),...
%         find(stru_grid(:,1)==stru_quad_up_plus(i,4)),...
%         find(stru_grid(:,1)==stru_quad_up_plus(i,5)),...
%         find(stru_grid(:,1)==stru_quad_up_plus(i,6)));
%     p=p+1;
% end
% for i=1:size(stru_quad_down_plus,1)
%     if stru_quad_down_plus(i,1)==-1
%         break
%     end
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_down_plus(i,3)),...
%         find(stru_grid(:,1)==stru_quad_down_plus(i,4)),...
%         find(stru_grid(:,1)==stru_quad_down_plus(i,5)),...
%         find(stru_grid(:,1)==stru_quad_down_plus(i,6)));
%     p=p+1;
% end
% fclose('all');
% !del 159.pch;
% !del 153.pch;
% stress_pT=stress;
% save result\stress_pT.mat stress_pT
%% 总应力（弃案）
% copyfile('result\103_end.pch','103.pch');
% WW=W+D_T(:,2:7);
% 
% f=fopen('spc2.bdf','w');
% for i=1:size(stru_grid,1)
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),1,WW(i,1));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),2,WW(i,2));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),3,WW(i,3));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),4,WW(i,4));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),5,WW(i,5));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),6,WW(i,6));
% end
% fclose(f);
% 
% eval(['!',Npath,' 101_post.bdf']);
% check_f06('101_post.f06');
% delete('101_post.f04');
% delete('101_post.f06');
% delete('101_post.log');
% delete('101_post.xdb');
% delete('101_post.DBALL');
% delete('101_post.MASTER');
% delete('101_post.IFPDAT');
% !del 103.pch;
% 
% f=fopen('101_post.pch');
% z=1;
% stress=zeros(0,2);
% while true
%     tline1=fgetl(f);
%     if length(strfind(tline1,'HEXA'))==1
%         while true
%             tline1=fgetl(f);
%             if tline1(1)=='$'
%                 break
%             end
%             stress(z,1)=sscanf(tline1(1:16),'%d');
%             for j=1:3
%                 tline1=fgetl(f);
%             end
%             stress(z,2)=sscanf(tline1(60:72),'%f');
%             z=z+1;
%             for k=1:60
%                 fgetl(f);
%             end
%         end
%         break
%     end
% end
% fclose(f);
% !del 101_post.pch;
% stress(:,2)=stress(:,2)/1000;
% f=fopen('plot\stress_pT.plt','w');
% fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
%     ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)+size(stru_quad_up_plus,1)+size(stru_quad_down_plus,1)-2),...
%     ', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
% l=1;
% for i=1:size(stru_grid,1)
%     fprintf(f,'%e ',stru_grid(i,2));
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% l=1;
% for i=1:size(stru_grid,1)
%     fprintf(f,'%e ',stru_grid(i,3));
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% l=1;
% for i=1:size(stru_grid,1)
%     fprintf(f,'%e ',stru_grid(i,4));
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% l=1;
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_up)
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_down)
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_up_plus(:,1))
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_down_plus(:,1))
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% p=0;
% for i=1:length(P_aero_up)
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),2)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),6)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),10)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),14)));
%     p=p+1;
% end
% for i=1:length(P_aero_down)
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),2)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),6)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),10)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),14)));
%     p=p+1;
% end
% for i=1:size(stru_quad_up_plus,1)
%     if stru_quad_up_plus(i,1)==-1
%         break
%     end
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_up_plus(i,3)),...
%         find(stru_grid(:,1)==stru_quad_up_plus(i,4)),...
%         find(stru_grid(:,1)==stru_quad_up_plus(i,5)),...
%         find(stru_grid(:,1)==stru_quad_up_plus(i,6)));
%     p=p+1;
% end
% for i=1:size(stru_quad_down_plus,1)
%     if stru_quad_down_plus(i,1)==-1
%         break
%     end
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_down_plus(i,3)),...
%         find(stru_grid(:,1)==stru_quad_down_plus(i,4)),...
%         find(stru_grid(:,1)==stru_quad_down_plus(i,5)),...
%         find(stru_grid(:,1)==stru_quad_down_plus(i,6)));
%     p=p+1;
% end
% fclose(f);
% stress_pT=stress;
% save result\stress_pT.mat stress_pT
%% 热致应力（弃案）
% copyfile('result\103_end.pch','103.pch');
% WW=D_T(:,2:7);
% 
% f=fopen('spc2.bdf','w');
% for i=1:size(stru_grid,1)
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),1,WW(i,1));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),2,WW(i,2));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),3,WW(i,3));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),4,WW(i,4));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),5,WW(i,5));
%     fprintf(f,'SPC     %8d%8d%8d%.4f\n',1,stru_grid(i,1),6,WW(i,6));
% end
% fclose(f);
% 
% eval(['!',Npath,' 101_post.bdf']);
% check_f06('101_post.f06');
% delete('101_post.f04');
% delete('101_post.f06');
% delete('101_post.log');
% delete('101_post.xdb');
% delete('101_post.DBALL');
% delete('101_post.MASTER');
% delete('101_post.IFPDAT');
% !del 103.pch;
% 
% f=fopen('101_post.pch');
% z=1;
% stress=zeros(0,2);
% while true
%     tline1=fgetl(f);
%     if length(strfind(tline1,'HEXA'))==1
%         while true
%             tline1=fgetl(f);
%             if tline1(1)=='$'
%                 break
%             end
%             stress(z,1)=sscanf(tline1(1:16),'%d');
%             for j=1:3
%                 tline1=fgetl(f);
%             end
%             stress(z,2)=sscanf(tline1(60:72),'%f');
%             z=z+1;
%             for k=1:60
%                 fgetl(f);
%             end
%         end
%         break
%     end
% end
% fclose(f);
% !del 101_post.pch;
% stress(:,2)=stress(:,2)/1000;
% f=fopen('plot\stress_T.plt','w');
% fprintf(f,['TITLE     = "Finite-Element Data"\nVARIABLES = "X"\n"Y"\n"Z"\n"Element UserID"\nZONE T="Step 0 Incr 0"\n STRANDID=1, SOLUTIONTIME=0\n' ...
%     ' Nodes=',num2str(size(stru_grid,1)),', Elements=',num2str(size(aerobox,1)+size(stru_quad_up_plus,1)+size(stru_quad_down_plus,1)-2),...
%     ', ZONETYPE=FEQuadrilateral\n DATAPACKING=BLOCK\n VARLOCATION=([4]=CELLCENTERED)\n DT=(DOUBLE DOUBLE DOUBLE DOUBLE )\n']);
% l=1;
% for i=1:size(stru_grid,1)
%     fprintf(f,'%e ',stru_grid(i,2));
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% l=1;
% for i=1:size(stru_grid,1)
%     fprintf(f,'%e ',stru_grid(i,3));
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% l=1;
% for i=1:size(stru_grid,1)
%     fprintf(f,'%e ',stru_grid(i,4));
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% l=1;
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_up)
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_down)
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_up_plus(:,1))
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% for i=1:size(stress,1)
%     if ismember(stress(i,1)-200000,stru_quad_down_plus(:,1))
%         fprintf(f,'%e ',stress(i,2));
%     else
%         continue
%     end
%     if l==5
%         fprintf(f,'\n');
%         l=1;
%     else
%         l=l+1;
%     end
% end
% fprintf(f,'\n');
% 
% p=0;
% for i=1:length(P_aero_up)
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),2)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),6)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),10)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_up(i),14)));
%     p=p+1;
% end
% for i=1:length(P_aero_down)
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),2)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),6)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),10)),...
%         find(stru_grid(:,1)==aerobox(aerobox(:,1)==stru_quad_down(i),14)));
%     p=p+1;
% end
% for i=1:size(stru_quad_up_plus,1)
%     if stru_quad_up_plus(i,1)==-1
%         break
%     end
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_up_plus(i,3)),...
%         find(stru_grid(:,1)==stru_quad_up_plus(i,4)),...
%         find(stru_grid(:,1)==stru_quad_up_plus(i,5)),...
%         find(stru_grid(:,1)==stru_quad_up_plus(i,6)));
%     p=p+1;
% end
% for i=1:size(stru_quad_down_plus,1)
%     if stru_quad_down_plus(i,1)==-1
%         break
%     end
%     fprintf(f,'%d %d %d %d\n',find(stru_grid(:,1)==stru_quad_down_plus(i,3)),...
%         find(stru_grid(:,1)==stru_quad_down_plus(i,4)),...
%         find(stru_grid(:,1)==stru_quad_down_plus(i,5)),...
%         find(stru_grid(:,1)==stru_quad_down_plus(i,6)));
%     p=p+1;
% end
% fclose(f);
% stress_T=stress;
% save result\stress_T.mat stress_T
%% Show Show Way
fclose('all');
end